clc;
clear all;
syms p
x = input('Input x values as an array : ');
y = input('Input y values as an array : ');
x0 = input('Enter the x for which y is to be found : ');
%Finding the length of x
n = length(x);
%Divided difference formula
for i = 1:n
 D(i,1) = y(i);
end
for i = 2:n
 for j = 2:i
 D(i,j)=(D(i,j-1)-D(i-1,j-1))/(x(i)-x(i-j+1));
 end
end
fx = @(p) D(n,n);
for i = n-1:-1:1
 fx = fx*(p-x(i)) + D(i,i);
end
fx = simplifyFraction(fx)
fprintf('Newtons iterated value using divided difference : f(%f) = %f \n',x0,
subs(fx,p,x0));